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Abstract. 

Fourier series of smooth, non-periodic functions on [—1,1] are known to exhibit the Gibbs phenomenon, and exhibit 
overall slow convergence. One way of overcoming these problems is by using a Fourier series on a larger domain, say 
[—T, T] with T > 1, a technique called Fourier extension or Fourier continuation. When constructed as the discrete least 
squares minimizer in equidistant points, the Fourier extension has been shown shown to converge geometrically in the 
truncation parameter N. A fast 0{N log^ N) algorithm has been described to compute Fourier extensions for the case 
where T = 2, compared to 0{N^) for solving the dense discrete least squares problem. We present two 0(7Vlog^ N) 
algorithms for the computation of these approximations for the case of general T, made possible by exploiting the 
connection between Fourier extensions and Prolate Spheroidal Wave theory. The first algorithm is based on the explicit 
computation of so-called periodic discrete prolate spheroidal sequences, while the second algorithm is purely algebraic 
and only implicitly based on the theory. 


1. Introduction. Fourier series are a good choice for the approximation of a smooth periodic 
function on a bounded interval. They offer exponential convergence, good frequency resolution, and 
the approximation can be computed numerically via the FFT. However, when the function is smooth 
but non-periodic, the exponential convergence of a Fourier series over the interval is lost, and ringing 
artefacts known as the Gibbs phenomenon are introduced. 

The Fourier extension technique (FE) [5, 6, 7, 8] aims to transfer the desirable properties of 
Fourier series for periodic functions to the non-periodic case. The principle is to approximate a non¬ 
periodic function that is defined on [—1,1] by a Fourier series that is periodic on [—T, T]. While the 
approximation may vary wildly in [—T, —1[ and JliT], under certain conditions it is guaranteed to 
converge exponentially to the original function within the interval. An illustration is shown in Figure 
1.1, where the extension is seen to agree closely with the given function on [—1,1]. Outside this interval, 
the extension is arbitrary, and in most cases defined by the solution method. 

The main difficulty with this technique is the ill-conditioning of the restricted Fourier basis. Nu¬ 
merically, this leads to ill-conditioned linear systems, which are difficult to solve efficiently. 



Fig. 1.1: A periodic extension to [—T,T] of a smooth function on [—1,1]. 


These constructions have been known in embedded or fictitious domain methods for general bases. 
A study on the approximation properties of extensions in the Fourier basis by Boyd [5] revealed that 
inside the smaller interval, the extension can be exponentially converging to the function. Further, 
he proposed the truncated singular value decomposition as a robust method for computing extensions 
from equispaced data. The resulting scheme, named FPIC-SU, can compute extensions of functions 
in the smaller interval, that are exact almost up to machine precision. The convergence rate is only 
limited by the smoothness of the function. Bruno et. al. used the same principle when approximating 
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surfaces by Fourier series on extended domains [7]. This provided a starting point for the very efhcient 
FC-Gram method [9, 25, 4]. 

Exponential error convergence of the FE problem was proven in [18] when inverting the Grammian 
matrix of the continuous least squares problem. Later, Adcock et. al extended the convergence analysis 
to the discrete least squares problem on equispaced data. At least superalgebraic convergence was 
proven for analytic functions, when using the truncated SVD. 

From an implementation point of view, the cost of the full SVD required for the FPIG-SU is 
prohibitively large. Recently, an FE algorithm was introduced by Lyon [23] that computes extensions 
in 0{Nlog'^ N) time. However, this algorithm only produces extensions of double length, i. e. it is 
unique to the case T = 2. Due to the reliance on symmetries only present when T is a power of 2, the 
algorithm cannot easily be extended to arbitrary T. 

In this paper we present fast algorithms for the computation of Fourier extensions of arbitrary 
extension length. One argument for varying the parameter T is found in the resolution power of 
the extension. The number of degrees of freedom per wavelength to represent an oscillatory function 
approaches the optimal value of 2 as T approaches 1 from above [1]. When T = 2, that number has 
already doubled. Other arguments may be performance related, as one may for example tune the 
length of the FFTs that are used in the computations, or a restriction on the data may be found in 
the application itself. 

Our algorithms stem from connecting the FE problem with classical results from signal processing 
theory. We state how it is essentially equivalent to the problem of bandlimited extrapolation. Gentral 
in this discussion are the so-called Prolate Spheroidal Wave functions, originally introduced at Bell 
labs in the 1970s [30, 21, 22]. The study of these functions and their special properties has been an 
active domain in signal processing since. 

The connection with Prolate Spheroidal Wave theory leads to explicit formulations for eigenvectors 
of the FE problem. Analysis of the EE problem learns that the relevant ill-conditioning can be captured 
using just (!l(loglV) of these vectors. Gombined with a fast solver for the remaining well-conditioned 
problem this leads to 0{N\og^ N) solvers. 

We explore two approaches in detail. In the first approach, a set of (!l(logV) discrete prolate 
spheroidal sequences is explicitly computed. This is based on their known properties: the vectors are 
known to be eigenvectors of a tridiagonal matrix. The second approach is purely algebraic and is not 
explicitly based on prolate spheroidal wave theory. Hence, this approach is more generally applicable. 
Both approaches have comparable performance, with a slight edge for the first approach. 

1.1. Overview of the paper. We formally state the FE problem in §2 and summarize relevant 
previous results. A concise overview of Prolate Spheroidal Wave functions is presented in §3, including 
subsequent work on discrete variants. The connection between these discrete variants and the EE 
problem is used to obtain fast algorithms in §4. Finally, we present numerical experiments to illustrate 
the numerical performance of these algorithms in §5. 

2. Fourier extensions. 

2.1. Problem formulation. For the remainder of this paper we will focus on infinitely differ¬ 
entiable non-periodic functions / on the interval [—1,1]. Other intervals are easily dealt with through 
affine transformations. 

The approximation is constructed on the extended interval [—T, T], with T the extension param¬ 
eter. For ease of notation denote 






Then 
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is the set oi N = 2n + 1 Fourier basis functions on the extended interval. The Fourier extension 
problem is now formalised as finding an approximation FAr(/) snch that 


FN{f) = min ||/-g||[-i,i]. ( 2 . 1 ) 

g^yN 

We note that the infinite set Goo is a basis for on the extended interval [—T, T], but it constitutes 
a frame on [—1,1] [12], i.e. it is redundant. This reflects the fact that a function can be extended in 
different ways. Once truncated to finite N, Gn is actually a basis on [— 1 , 1 ], albeit a very ill-conditioned 
one. We refer to [2] for a detailed discussion of this point of view. 

In an implementation context, it is more natural to consider the search for the Fourier coefficients 
oG C^, 

n 

a= min II/- V Ck<pk\\[-i,i]- 

k——n 

Depending on the norm to be minimised, the extension problem takes on a different formulation: 

Discrete Fourier extensions. Most practical applications provide information about / as samples 
in a predefined set of points. The norm in (2.1) is then conveniently replaced by a discrete least squares 
norm 

pNif) = - g{xi)f , xi G [-1,1]. (2.2) 

seSjv 

Previous work distinguishes between uniform sampling, and an optimal sampling set called mapped 
symmetric Chebychev nodes [1]. The focus of this paper is on the case where the sampling points are 
uniform, 

I , 

Xl = —, I = —TO, . . . , TO, 

TO 

for a total of M = 2m + 1 points. Typically, some oversampling 7 is used so M = jN > N. Including 
appropriate normalisation, this leads to the matrix formulation that is used throughout the rest of the 
paper. Let 

1 .^rkl I 

Ai,k = = -^==e'rm , bi = —=f{xi). (2.3) 

V2Tm Vm 

Note that matrix-vector products involving A can be performed very efficiently using FFTs of length 
L = 2Tm. Therefore L is assumed to be integer for the remainder of this paper. 

The Fourier coefficients are then found by collocation, through solving the rectangular system 


in a least squares sense. Note that due to the nature of the problem we are not interested in a specific 
solution a, just one of the possibly many solutions that lead to a small residual ||^a — 6 ||. Solving this 
problem efficiently is the focal point of this paper. 

Continuous Fourier extensions. When the norm is the regular /!^-norm over [—1,1], 

Fnif) = min ||/- ffibj-i.i] 

g^yN 

is known as the continuous Fourier extension. The resnlting least-squares problem can be solved by 
formulating the Grammian matrix. Let 


Ai,j — 


4>i(x)4>j (x)dx = 


K* - 3 ) 


bi= f(l>i(x)dx. 


Then the Fourier coefficients follow from 


Aa = b. 
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2.2. Convergence and stability. The usefulness of approximation schemes hinges on two prop¬ 
erties: the speed of convergence to the given function, and the stability of the required computations. 
Considerable effort has been put into quantifying these properties both analytically and numerically 
[18] [1] [2] [24]. Without going into too much detail, we recap the most important results. 

First of all, a distinction should be made between the exact discrete and continuous FE solutions 
FN{f) and Fn^J), and their computer-implemented counterparts. Computing the exact solutions 
is known to be unstable, as they can grow unbounded outside the interval of interest. Numerical 
algorithms however will never compute these exact solutions. Due to regularisation, the numerical 
FEs GAr(/) and GAr(/) are more stable, while maintaining the desired convergence behaviour. In [2], 
the truncated SVD was used as a model scheme for solving the ill-conditioned systems (2.3) and (2.5). 
Given the Singular Value Decomposition of the FE matrix A = USV, the solution to Aa = b is 


= VS'^U'b, 



Si,i > T 

otherwise. 


(2.7) 


The truncation parameter r then acts as a regularisation parameter. It is usually chosen close to the 
machine precision. 

Stability. Following [2], stability is defined in terms of the absolute condition number of the FE 
mapping 


«(Fv) = sup{||E^( 6 )|| : 6 G ll&ll = 1}, 

where, with slight abuse of notation, Jbv( 6 ) is the solution to the FE problem with right hand side b, 
and jj - II is the regular P norm over [—1,1]. This condition number can be computed for the continuous 
and discrete FEs, both the exact and numerical versions. It was shown to grow exponentially in N for 
the exact solution to the continuous and discrete FE problem. 

When looking at the numerical FEs the situation changes considerably. The condition number 
of the numerical continuous FE mapping is k{Gm) ^ l/v^- The condition number of the numerical 
discrete FE k{Gm) is dependent on a constant 0 < a{j]T) < 1, independent of N, that satisfies 
a{--f;T) —>■ 0 as 7 —>■ oo for fixed T. It is given by k{Gn) ^ \fN G N. This means that for 

a sufficiently large oversampling factor 7 , the condition number of the numerical FE mapping can be 
made reasonably close to 1 . 

Meanwhile, the condition number of the matrices A and A grows exponentially as V —>■ 00 , where 
M > N for the discrete FE. This is surprising, given the good condition of the FE mapping. It can be 
understood by noting that extensions with small coefficient norm and small residual are guaranteed 
to exist. The numerical algorithms will steer clear of the unstable exact solution, and instead return 
one of these alternatives. For a full exposition on the stability of FE calculations, see [2]. 

Convergence. Concerning convergence, results in [2] are valuable only for functions that are ana¬ 
lytic in a region I 2 (p*) of the complex plane and continuous on its border. This region is a Bernstein 
ellipse under a transformation that allows Fourier extensions to be understood as polynomial approx¬ 
imations. For such functions /, the exact continuous and discrete FEs converge geometrically, with a 
speed 


||/-Fv(/)|| <c/p-^. 

Here p = min{p*,E(r)} and cy is proportional to max3,gx)(p) l/(a^)l- F{T) is known as the Fourier 
extension constant and is given by cot^ (^)- Note that for the exact discrete FE, there is an added 
requirement of scaling M as to avoid the Runge phenomenon. 

Under the same analyticity conditions, the error decay of the numerical counterparts Gjq and Gn 
can be broken down into several subregions: 

1 . li N < N2, where N2 is a function-independent breakpoint, ||/ —GAr(/)|| converges or diverges 
exponentially fast at the same rate as the exact solution. 
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2. When N < Nq (continuous) or N 2 < N < Ni := 2 Nq (discrete), where Ng is another function- 
independent breakpoint depending, both ||/ — GAr(/)|| and ||/ — G'Ar(/)|| decay like p~^■ 

3. When N = Nq or N = Ni, the errors are approximately 

\\f-GNoif )\\^ 11 / - Gn, (/)|| « 

where c/ is as before, and df = ^ 

4. When N > Nq or N > Ni, the errors decay at least super algebraically fast down to maximal 

achievable accuracies of order y'r and respectively. 

This behaviour of the error offers insight into the usability of Fourier extensions. A first observation 
is that the continuous FE is limited to a maximal achievable accuracy of ^/r. Coupled with the need 
to compute Fourier integrals to compose the right hand side b in (2.5), this makes the continuous FE 
unfit for practical use. However, the algorithms presented in §4 for the discrete EE can be adapted to 
this context with little extra effort. This is documented in §4.3. 

On the contrary, the numerical discrete EE guarantees convergence up to a certain power of r. By 
varying this cutoff, the oversampling and the extension length T this maximum achievable accuracy 
can be made very close to the machine precision. 

2.3. Influence of the extension length. The main contributions of this paper are algorithms 
that add flexibility in the choice of extension length T. The increased resolution power was already 
cited as an argument to reduce T, but it is important to be aware of the possible consequences. 
Therefore, in this section we summarize the influence of this parameter on convergence, resolution 
power, and conditioning of the EE problem. 

First note that the Fourier extension constant E{T) grows with T. For functions analytic in a 
sufficiently large region, the convergence rate p is limited by this constant. Increasing T thus increases 
the convergence rate, and vice versa. 

The resolution power of a scheme, first studied by Gottlieb and Orszag [15] is a measure of the 
amount of point samples needed to resolve an oscillatory function to a certain precision. Let 

n{oj,5) = mm{NeN : - Fjv(e*’^‘^)||oo < 5}, w > 0, 

for some small S. Then Ejy has a resolution constant r if 

7Z(uj, d) ~ ruj, u) —>■ 00 . 

For regular Fourier series, this constant has the optimal value 2. 

A theoretical argument shows that for the continuous FE this resolution constant increases with 
T [1]. More specifically. 


r(T) < 2Tsin (^), Tg(1,oo). 

Thus, for T « 1 the resolution constant r{t) « 2T is close to optimal. When T tends to infinity, 
r(T) ~ TT. It is even possible to optimally balance convergence speed with resolution power when 
aiming for a predetermined accuracy etoi- This is achieved by varying T with N , specifically 

T(fV,eto/) = ^ (^arctan(etoi)5^^ . (2.8) 

Although no equivalent analysis exists for the discrete EE, there have been several attempts to 
determine EE parameters that are in some sense optimal. In [ 8 ] Bruno suggests the values T = 2 and 
7 = 2 as a general rule of thumb, but at the same time notes that the optimal parameters are heavily 
function dependent. Note that increasing both the extension length T and the oversampling 7 will 
likely increase the resolution constant r. Especially since it was shown in [1] that the limit r(T) ^ tt 
as T —>■ 00 no longer holds for a discretised EE, instead the resolution constant grows as r(T) ~ 2T. 
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Even though this was only observed for data points distributed as a variant of Chebychev points, it is 
indicative that the resolution constant for T = 2,7 = 2 will be considerably above the optimal value. 

The precise interplay between T and 7 on the one hand, and resolution power and conditioning on 
the other hand was studied in detail by Adcock and Rua [3]. They found that the condition number 
k{Gn) of the equispaced discrete FE depends only on the product T 7 . Increasing either will lower 
the condition number. Thus as long as 7 is increased or decreased accordingly when varying T, the 
conditioning of the FE mapping remains constant. This is cited as an argument to limit T to 2, to 
profit from the at the time only available fast FE algorithm. 

Furthermore, the resolution constant is also dependent on the product T 7 , growing as r(T) ~ T 7 . 
This illustrates the tradeoff between resolution power and conditioning. Numerical experiments in [3] 
showed that by allowing the condition number to grow from k « 10 to k « 100, the resolution constant 
was halved, while further increasing n had very little additional value. However, it should be noted 
that these experiments were only carried out for T = 2. Lifting the restriction on T may thus offer 
more flexibility in finding a balance between resolution power and conditioning. 

An interesting open problem raised in [3] is the possibility to vary T with M, to achieve optimal 
resolution power in a manner similar to (2.8). Due to the lack of a fast algorithm, any gains from 
varying T were considered of limited practical usability compared to the fast algorithm for T = 2. The 
fast algorithms presented in this paper warrant a closer look at the possible benefits from this method. 

3. Prolate spheroidal wave functions and discrete variants. A long-standing problem in 
signal processing theory is that of bandlimited extrapolation. The problem is, assuming some portion 
of a bandlimited signal is known, to accurately predict the missing data. In a first subsection we explain 
how the FE problem is a specific variant of this problem. The subsequent sections then explore the 
theory of Prolate Spheroidal Wave functions, that plays a major role in bandlimited extrapolation, for 
further use in the FE algorithms. 

3.1. Discrete Fourier extensions and bandlimited extrapolation. The discrete Fourier 
extension is closely related to discrete band limited extrapolation, i.e., to the problem of reconstructing 
a discrete bandlimited signal from a number of data samples. Simply put, we are looking for a vector 
y with discrete Fourier transform Y such that 

y[k] ^ f[k] k G SM,t 
Y[l] =0 Sn,u., 

where / is the sample data, and SM,t and Sn,ui are sampling sets in the discrete time and frequency 
domains, of sizes M and N respectively. Depending on the problem parameters the solution might not 
be unique. In this case an additional minimal solution norm constraint can be added [19]. 

A popular method is the Papoulis-Gerschberg algorithm [14] [26]. It uses a two-step iteration 
process to alternate matching the given data and complying with the frequency constraints. Variants 
that use the conjugate gradients and related methods to speed up the iteration process tend to perform 
reasonably well numerically [31]. They operate at a cost of 0{N log N) operations per iteration, where 
the number of iterations scales with the bandwidth of the signal. 

Besides these iterative methods, considerable attention was given to direct methods for solving 
this discrete problem. One such method is the one proposed by Jain and Ranganath [19], where both 
the time and frequency sampling sets are contiguous, i. e. St = —m,... ,m and S^i = —u,...,n. They 
commence by writing the data as a function of the unknown coefficients: 


f = Jy = DmBnV 

5i,j ) 




D 1 ■ (P “ 

Bn,p,i = j 2 ^ exp 1 --- 


k— — n 



(3.1) 

(3.2) 

(3.3) 
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Here, is a L x L circulant matrix that represents a discrete low-pass filter, and Dm is a M x L 
selection operator. Sij is the kronecker delta. D'j^ is an extension operator that pads a sequence of 
length M with zeros to length L. 

The direct methods of Jain and Ranganath then consist of solving Jy = f in a least squares sense 
by formulating the normal equations 


J'Jy = J'f- 


Since J' J is symmetric positive definite and Toeplitz, the Levinson-Trench algorithm can be applied 
to compute the inverse of J'J in 0{L'^) operations. They also suggested another approach, using the 
singular value decomposition of J to solve the least-squares problem. The resulting singular vectors 
were named periodic discrete prolate spheroidal sequences (P-DPSSs), after the prolate spheroidal 
wave functions (PSWFs) arising in continuous bandlimited extrapolation. 

A closer look at the matrix A from the discrete Fourier extension (2.3) makes the relation with 
(3.1) apparent. Adopting the notation for the DFT length L = 2Tm, 

, . 1 ■{P~ 9)27rA: 

[AA )pq = 2^ exp 1 - - -, p,q = -TO, ..., TO 

k——n 

= {DMBND'M)pq = {J j')pq, (3.4) 

where the last line follows from the idempotency of Bm- Consequently, A and J share the same left 
singular vectors, and the same singular values. Essentially, discrete Fourier extension is a reformulation 
of the bandlimited extrapolation problem with the N frequency coefficients as unknowns, instead of 
the extrapolated signal. The focus has also shifted from determining an accurate extension to the 
approximation of the given data samples. 

The next sections concern the PSWFs and P-DPSSs, and how they are natural solutions for the 
bandlimited extrapolation problem. The groundwork for this theory was detailed in a series of papers 
by Slepian, Landau and Poliak from 1961 onwards [30, 21, 22, 27, 28]. An overview is given in [29]. 

3.2. Prolate Spheroidal Wave Functions. Denote by f{x) and a function in C? and its 
Fourier transform, so that 

/ OO pOO 

-OO J —OO 

The time- and bandlimiting operators D and B are then defined as 

Vf{x) = fix) = I I""! J I Bf{x) = r (3.5) 

(0 \x\ > T J_Q 


which project onto and PWq, the Paley-Wiener space of bandlimited functions, respectively. 

Note that the bandlimiting operator can also be written as 


BfX) 


sin(2^n(x-.)) ^^ 

, Tr{x — s) 


The Heisenberg-Gabor limit states that no function can be simultaneously concentrated in both 
time and frequency, and so Vf : jjHU/jj < U/]]. However, one can look for nearly-invariant functions 
under this operator, functions for which j jHIJ/J]/]]/]] is as close to 1 as possible. 

These are the eigenfunctions of the operator BV, i. e. the solutions of the integral equation 

> , / N ,/ Nsin(27rH(a; - s)) 

X'lpix) = / 7/>(s)-— ^^ds. (3.6) 

J-T n{x-s) 
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Slepian and collaborators showed that this equation is solvable only for select values of A, a 
countably infinite set 1 > Aq > Ai > • • • > 0. The corresponding eigenfunctions tpi were named Prolate 
Spheroidal Wave functions. The naming stems from the curious observation that these functions are 
solutions to the spheroidal wave equation 

“ 1^) ^ - (27rf}T)^ (3.7) 

This is a Sturm-Liouville equation with a set of unique eigenvalues ••• > 6i-i > 9i > 0i+i > ... 
corresponding to the functions 

Since B and T) are idempotent operators, it is convenient to consider the "(pi eigenfunctions of the 
Hermitian operator BT^B. Timelimiting both sides of (3.6), the timelimited functions tpi = 'D'lpi are the 
eigenfunctions of the Hermitian operator 'DB'D, with corresponding eigenvalues Ai. The term Prolate 
Spheroidal Wave function is used for both the %pi and the ipi. 

As eigenfunctions of a Hermitian operator, the tpi and ipi are orthogonal 

ipi{x)ipj{x)dx = 5ij J 'tpi{x)'tpj{x)dx = XiSij, 

and they are complete in PWq and respectively. The Prolate Spheroidal Wave functions 

thus form an orthogonal base for PWq, and the eigenvalue Ai represents the fraction of energy of ipi 
contained in [—T, T]. 

This leads to a straightforward approach to continuous bandlimited extrapolation. Let / be a 
function segment in ■ Then 



9 = 




(3.8) 


is a bandlimited function that agrees with / in the interval due to the completeness of the ipi in • 

Furthermore, when truncating the sum, the first terms have the largest eigenvalues and capture the 
most energy inside the interval. 

Using both (3.6) and (3.7), the PSWFs were further shown to have the following properties [30]: 
i The Ipi are eigenfunctions of the finite Fourier transform. 


P‘^'^*^'tpn{t)dt = : 



ii The eigenvalues Ai cluster near 1 for low values of i, and decay exponentially after a set breakpoint 


Ai « 1, i <C 4HT and Ai « 0, i ^ 4flT. 

The width of the region where Ai G (e, 1 — e) grows as log(HT). 

The first property illustrates the symmetry of time and frequency domain in this problem. The 
second property shows that a truncated approximation (3.8) will require « 4flT terms. The increasing 
irregularity of the (pi coupled with the exponential decay of the Ai ensures that g will converge rapidly 
after this breakpoint for a sufficiently regular /. 

3.3. Periodic Discrete Prolate Spheroidal Sequences. There are several possible discreti¬ 
sations for PSWFs. The most well-known is the one proposed by Slepian, where the Fourier transform 
is replaced with the discrete-time Fourier transform 

G{f) = E 9[n]e-^^^f\ 9[n] = - / G{f)P^-f-df. 
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The resulting generalisations are frequency domain functions on [—1/2,1/2] commonly referred to as 
discrete prolate spheroidal wave functions (DPSWFs), and infinite sequences in the time domain called 
discrete prolate spheroidal sequences (DPSSs). The properties from section 3.2 carry over, apart from 
some loss of symmetry between time and frequency domains [28]. 

Another approach was proposed independently by Jain and Ranganath [19] and Griinbaum [17], 
and elaborated on in [33]. They replace the Fourier transform with the DFT for sequences of length 
L: 


L-l L-1 

Hk = y\ h[n] = y V 

ra =0 fc =0 

The discrete analogues of the time- and bandlimiting operators are given by the matrices (3.2) and (3.3) 
from the discrete band limited extrapolation problem. Similar to the continuous PSWFs, the periodic 
discrete prolate spheroidal sequences (P-DPSS) (jii are defined as the eigenvectors of the hermitian 
matrix operator BjqD'j^DMBN, and their limited versions (pi as the eigenvectors of DmBnDm 

BnD'j^DmBnPi = Vipi DMBND'j^pi = Vipi. (3.9) 


Since DmBnD'j^ is not of full rank when M > N, an additional requirement that Vi ^ 0 is imposed, 
following [33]. Consequently, there are only {min(Af, M)} P-DPSSs. 

The pi properties are similar to those of PSWFs: 

i If M > A^, the pi are complete in the space of bandlimited sequences spanned by Rat. 

ii If M < N, the pi are complete in 

iii The pi are doubly orthogonal 

Pi * Pj — ^ij( Bpi * Dpj — ViSij 

iv The P-DPSSs are eigenvectors of the finite DFT, but with the roles of M and N interchanged 

D]\[FpM,N,i = pN,M,i- 


Here, F is the DFT matrix of size L. 

V Among sequences of length L, with frequency support in [0,7Ti.] U [L — m, L], po is the most con¬ 
centrated in [0,M — 1]. Among sequences of equally limited frequency support orthogonal to po, 
pi is the most concentrated in [0,M — I], and so on. 
vi Like the eigenvalues of the PSWFs, the eigenvalues Vi are distinct and cluster exponentially near 
1 and 0, in that 


Vi 


. . nm 

I, i <C —-— and Vi 


0,1 > 


NM 

L 


The width of the plunge region where Vi G (e, 1 — e) grows as \og{NM/L) [32]. 
vii The P-DPSSs satisfy a second order difference equation [33] 

bk(i>i[k - 1] + Ck(pi[k] + bk+i(pi[k + 1] = 6i(pi[k], k = 0,...,M-l, (3.10) 


with coefficients 


bk 





— COS 



\ . /7r(M — k)\ 

(2k+l-M)\ 

L ) 
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Equation (3.9) clarifies their importance in both the discrete band limited extrapolation problem, 
and discrete Fourier extensions. From (3.4) it is obvious that the eigenvectors of AA' and JJ' are the 
limited P-DPSSs (^i. Regarding the singular value decompositions of A and J, this determines the left 
singular vectors as (/>i and the associated singular values as For J, following 

J' J = = BnDmBn, 

the right singular vectors are seen to be full P-DPSSs (/)i. Using these relations, Jains SVD method 
results in an extrapolated sequence 


N-l ^ 

”] = E • f) 

i=0 

similar to the continuous extrapolation (3.8). 

The discrete Fourier extension matrix A differs from J in the right singular vectors. Similar to 
(3.4) 


A'A 



expi- 


. (p — q)27rk 


= DnBmDn- 


p,q= -n,...,n 


The right singular vectors are then again limited P-DPSSs 4>i, where the parameters N and M have 
been interchanged. Using $ to represent the set of P-DPSSs and T for the diagonal matrix of eigen¬ 
values, the SVD of A is given by 

A = 

Calculating the coefficients of the discrete FE by truncating the SVD now corresponds to 


a = 




where ic is determined by the cutoff parameter. 


> T > 

This formulation of the FE together with the P-DPSS properties listed above leads to the fast algo¬ 
rithms in the next section. 

4. Algorithms. In this section we present two approaches to computing the truncated SVD 
solution to the discrete Fourier extension problem Ax = b. Both approaches are based on a specific 
division in easier subproblems, which is explained below. They differ only in how they solve one of 
these subproblems. Detailed algorithms follow in sections 4.1 and 4.2. 

From the previous section, the SVD of the FE matrix A is 

A = USV' = 

The cost of this full SVD is prohibitively large, growing as O(N^). To reach the performance goal 
of 0(N log N) operations, two subproblems are identified and solved in succession. The key to this 
subdivision is in the distribution of the singular values S, as shown in figure 4.1. Three distinct regions 
are visible. Following property (vi) of the P-DPSS, they are for some small cutoff r: 

• A region /„ := {cr : 1 > ct > 1 — r} where all singular values are 1 up to a tolerance r. This 
region contains approximately 0{NM/L) singular values. 
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Fig. 4.1: The subdivision of the spectrum of A into three distinct intervals, with cutoff parameter 
r = le — 14. Due to rounding errors, the eigenvalues in region Ij don’t decay past machine precision. 


• A region := {a : 1 — t > a > t} , also referred to as the “plunge region” in a more general 
context regarding truncated frames. The number of singular values in this region grows as 
C>(log {NM/L)). 

• A region := {a : t > a > Q} where the singular values further decay exponentially. 

Let the subdivision of the singular values and associated vectors be denoted by 


A 


[u^ Up U^] 




S/3 




The truncated SVD solution to the problem Ax = b with truncation cutoff at r is then 


X=Xa+Xp = VaT.^^U^ba + Vpl^^^Upbp, (4.1) 

where the inverse operator applies to the diagonal elements. The right hand side is split along the 
orthogonal spans of Ua, Up and U-y, i.e., 


= UM, 

and so on. Note that this solution method implicitly assumes bj to be below the required solution 
accuracy. If it is not, a large solution term x.y is needed, with ||a:-,.|| > yll&^H. In this case, we say the 
Fourier extension has not yet converged, and N should be increased. 

Equation (4.1) splits the problem into two orthogonal subproblems. The isolation of Vp, Sp and Up 
from the plunge region and efficient calculation of is where the two approaches given in sections 4.1 
and 4.2 differ. 

The subsequent calculation of Xa is then straightforward, based on the following observation: 

A'{b - bp) = + V^Ys^U’^b^ 

= V^Y-'^U'^b^ + 0{T) 

= Xa + 0{t) 

The b.y term, which is already assumed to be negligible, is fully eliminated by the additional 0{t) 
factor E..^. Noting that Ea = E“^ + 0{t) then yields x^ at the cost of a single fast multiplication with 
A' {0{NlogN)). 
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4.1. Approach 1: explicit projection onto prolate sequences. Due to the intrinsic connec¬ 
tion of the FE problem with DPSSs, it is possible to explicitly compute and V^. The second 

order difference equation from (3.10) implies that 


TM,N4'M,N,i = 


Tm,N 


Co bi 
bi Cl 


bN-i 
bN-l Cn-i 


(4.2) 


The P-DPSSs can thus be found as eigenvectors of a tridiagonal matrix. The that make 

up U are eigenvectors of Tm,n- The dual matrix Tn,m yields the right singular vectors V. With the 
P-DPSSs known, the original singular value is found in 0{N log N) operations as 


This approach is already in use for the regular DPSSs [16]. 

In this stage of the algorithm however, we are only interested in a subset of the singular values 
and vectors. Since the number of singular values in the plunge region is logarithmically small, we 
look to algorithms that calculate k eigenvalues and eigenvectors of a tridiagonal matrix T in 0{kN) 
operations [11]. Algorithms for this computation require as input for the desired eigenvalues 9i of T 
either: 

• A range [Ci , C 2 ]'■ the algorithm then finds all 6 : Ci < 9 < C 2 - 

• A index set {imin) • ■ • > *max}: the algorithm then finds all 9 : 9i^^^ < 9 < from the 

ordered set 0o > ^1 > • • • ■ 

The algorithms thus require knowing the 9j, or alternatively the indices j, that correspond to 
(Ti G Ijs- Denote the mappings between ai and 9j, and between their indices i and j, by 


9j = j = 




The description of GN,MiIp) is difficult, as very little is known of this mapping. The only known trait 
is monotonicity, a trait shared with the PSWF and DPSWF equivalents. This is already very helpful, 
since it implies that G'pf m(*) = *• 

Theorem 1. For M > N, the mapping Q'j^ is monotone, 

Vii,f2 : *1 > *2 ^ ^Af,M(*2)- 


Proof. The proof follows a mechanism used by Slepian in both [30, p. 61] and [28, §4.1]. The 
continnity of eigenvalues and eigenvectors as a function of a parameter, combined with a known ordering 
result for a specific value of this parameter, extends the known result to all parameter values. 

First we recall a similar result from the continuous Fourier Extension. Let A be the matrix of the 
continuous EE (2.5) with eigenvalues ai, and T the tridiagonaln matrix (4.3) with eigenvalues 9i 

>Co bi 
rf, &1 Cl 


bN-l 


bN-i 

CN-1 


Ci = 


N -I 


— i 


cos 


T’ 


h = 


i{N-i) 


(4.3) 


12 



Then the mapping Q' 


3=Q\i)^ 


= Sj 

^iA^i = CTi 


was proven to be monotone in [30, §4.1]. This result can be extended to the discrete FE case, since 
the limits of the entries of A' A and Tm,m for large m can be written in terms of the corresponding 
matrices of the continuous problem: 


lim (A'A)ij = 


sm 




Tr{i - j) 


^ = A: 


TT^fcliV — k) 27r^ - 
hm bk = - 

m—¥oo Lj^ 

lim Ck = cos^i^^( -k\ - 1 

m^oo 1 \ \ 2 


27T^ _ TT 

= —TrCk — cos — 
T 


27r^ - TT 

hm Tn,m = -pr^ ~ cos —I 

m^oo ± 


The last line ensures that the eigenvalues of Tjv^m are those of T under a linear mapping. Since this 
mapping preserves the ordering of eigenvalues, the theorem holds in the limit m —>■ oo. Further, when 
limited to the regime T > 1, 2m + 1 > A^, we have the following: 

1. The tridiagonal matrix TM,Ar(m) with diagonal elements Ck and subdiagonal elements bk 




sm 


Ck = — cos 


7r(fc — m) 

Tm 


7r(2m + 1 — A:) 
2Tm 
f ttN \ 
\2TmJ 


commutes with A'A for integer values of m. However, by a substitution as in [10, Thm. 4.7], 
it is easy to see that this relationship holds for any real m > n. 

2. A classical result states that the eigenvalues of a matrix are continuous as the matrix entries 
vary continuously in the parameter. Thus all 0i(m) are continuous. In general they may 
coincide with each other. However, since all subdiagonal entries are always non-zero, 

VI < A: < iV - 1 : &^ > 0, 


7M.Ar(m) is a so-called normal Jacobi matrix and such matrices are known to have distinct 
eigenvalues [13, Ch. 2.1]. As a result of this distinctness, the eigenvectors can be chosen to be 
continuous in m as well [20, Ch. 2 §5.3]. 

3. To prove similar statements for the matrix A'A, we note that the inductive proof in [33, 
Prop. 5] for the distinctness of eigenvalues of A'A is dependent only on it being symmetric, 
and the commuting Tm,n{'<ti) being a normal Jacobi matrix. The eigenvalues Xi(m) are thus 
continuous and distinct Vm > n. Then the eigenvectors can also be chosen continuous in m. 
Combining these statements, the distinctness preserves the relative ordering of the continuous eigenval¬ 
ues. The continuity of the eigenvectors relates the eigenvalues of A and Tjv^m through the mapping Q'. 
Since this mapping is monotone in the limit m —>■ oo, the continuity ensures the mapping is monotone 
for all m > n. D 

The required index set for the eigenvalues of Tm,n is exactly the index set corresponding to ai 
from the plunge region. From (vi) these are known to be centered around and their number 

grows as (!l(log A^). All that remains is to determine the constants Crain and Cmax so that 


CTj G 1/3 j G 



Crain log N, N 
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NM 

L 


Ca 




The minimum required index set with cutoff r = le — 16 for increasing N is shown in Figure 4.2, 
with the value N — as a dashed line. Experimentally, the choices Cmin > 6 and Cmax > 3 seem 
sufficient for all r > Cmach- The Vp can be obtained by refining A'(j)M,N,i as eigenvectors of Tm,m- 
Using a fast tridiagonal eigenvector algorithm, {7/3, S /3 and Vp can be computed in 0{NlogN) 
operations. The solution term 


xp = 


is then found in 0{N log^ N) operations. 

Remark 1. Monotonicity of the map Qn,m is also observed to hold for integer M smaller than 
integer N. Specifically, a variant of Theorem 1 holds for the M nonzero singular values ai. The same 
index set can thus be used to compute both Up and Vp with a fast tridiagonal eigenvalue algorithm. 
However, this altered algorithm is without proof. 



Fig. 4.2: The behaviour of the index set of the plunge region. The minimal and maximal index of the 
plunge region are shown as solid lines, for different values of T. The point NM/L, which is known to 
lie in the interval, is shown as a dashed line. 


Combining xp with the calculation of Xa 
A bare-bones version is given in Algorithm 1. 


described above leads to a fast (!l(A^log^ N) algorithm. 


Algorithm 1 Explicit projection on DPSSs 

U/3, Vp =Gig(T3V, Tjvf, {jTnin j Jmax }) 

S/3 = {U'pAVp) > T 
xp = Vp^-p^U'pb 
Xa = A'{b — Axp) 

X = Xa + Xp 


4.2. Approach 2: an implicit projection method. The second approach to calculating xp is 
more universal, or more generally applicable, since it depends solely on the steep singular value profile 
discussed earlier and illustrated in Fig. 4.1. As such, it is extensible to any ill-conditioned basis that 
exhibits a similar profile. 

The approach is based on the observation that multiplying both the FE matrix A and right hand 
side by a factor 


P = {AA' - I) 
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(4.4) 








isolates the problem to the plunge region. This is most easily seen from the SVD. With A = USV' we 
have 


P = - T)U', 


(4.5) 


and 


PA = [/(s^ - j:)v'. 


Note here that the mapping 


P{<t) = — a 

isolates the singular values from the plunge region since Vcr G {la U Ij} '-Pier) = 0(t). This way, PA 
preserves the singular vectors of just the plunge region, but with mapped singular values 

PA = U0m-j:p)v;, + o{T). (4.6) 

In theory, P is a square full rank matrix, and solving 

PAx = Pb (4.7) 

is equivalent to solving Ax = b. In practice, PA has a large numerical nullspace. Hence, PA is 
approximately low rank, with the rank increasing with the size of the plunge region. The combination 
of this low rank with a fast matrix-vector product allows random matrix algorithms to solve (4.7) very 
efficiently. 

Assume W a uniform random matrix of dimensions N x R, where R = C log N -\-D is a conservative 
estimate for the rank of PA. From the previous section, C = C'min + C'max > 9 is sufficient, with P ~ 10 
minimizing the chance of failure of the random matrix algorithm. The column space of PAW then 
approximates the column space of PA very well. 

Therefore, solving the following small linear system 

PAWy = Pb 


and letting 


xw = Wy 

one obtains a solution to (4.7) a cost of 0{MR^). It follows from (4.5) and (4.6) that xp is recov¬ 
ered exactly. On the other hand, this solution process introduces additional solution terms from the 
nullspace of PA. Write xw as xp +ra + r^. Then as before we calculate 

A'{b — Axw) = A'{ba — Axa — Ar^) + 0{t) 

= Xa-Ta- + 0 {t) 

= Xa-ra+ 0(t). 

Then x = xw + A'(b — Axw) boils down to 


X = Xa + Xp + r.y + 0{t). 

so that 11 Ax — 6|| = 0(t}. The total cost of this algorithm is again 0{N log^ N) operations. A step 
by step version is given in Algorithm 2. 
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Algorithm 2 Implicit projection on DPSSs 
PAWy = Pb 
xp = Wy 
Xa = A'{b — Axp) 

X = Xa + xp 


4.3. Adaptation for the continuous FE. As mentioned in section 3.3, the continuous FE also 
has numerous connections with Prolate Spheroidal Wave theory. The solution is most easily expressed 
in terms of the eigenvectors of A, the DPSSs 


A = 


The eigenvalues Ai have a similar profile to that of Figure 4.1. Our second approach thus applies 
immediately, and Algorithm 2 is well suited to solve the continuous FE problem. Note however that 
this does not eliminate the theoretical 0(^/T) error bound. 

Furthermore, these DPSS also satisfy a second order difference equation, different from (3.10). In 
particular, the matrix A commutes with the tridiagonal matrix (4.3). 

It follows that, with minor modifications. Algorithm 1 can also be used to solve the continuous 
FE problem. 

5. Numerical Results. In this section we apply the algorithms from the previous section to a 
number of test problems. All tests were performed in matlab , single threaded. The required fast 
matrix-vector products Ax and A'y were implemented using ffts, and the lapack routine dstevx was 
used for the tridiagonal eigenvalue problem. 

To show the validity of our algorithms for different values of T, experiments are carried out for 

Ti = 1.1, T2 = 2, Ta = 3.8. 

Following [3], the product L = 2Tm is held constant when varying T in order to maintain a fixed 
condition number. Thus, in our computations we use 

Ml « —M2, M3 « ^-^2, 


with the restriction that M is odd. 

5.1. Computational complexity. Figure 5.1a shows execution time for increasing degrees of 
freedom of the algorithm for different values of T. The figure confirms the theorized 0{N \og^ N) 
asymptotic complexity of our algorithms. It also shows execution speed is on par with the current 
fast algorithm by Lyon. Also, as in Lyons algorithm, the majority of the work is in computing a 
low-rank matrix decomposition related to A, in this case the middle part of the SVD. This can provide 
a significant speedup when multiple approximations are needed with the same parameters. 

Figure 5.1b compares our two approaches in more detail. The figure shows the execution time per 
degree of freedom. The difference between algorithms is very much dependent on implementation, but 
it highlights a trend seen for different values of T. Both algorithms appear to be slightly faster for 
larger T, which is somewhat counterintuitive. After all, the cost of the FFT operations is 0{L), which 
is held constant for all experiments. However, both approaches solve a subproblem with a cost that 
grows as 0{M). This is either the tridiagonal eigenvalue problem with a matrix of size M x M, or a 
linear solve of an M x logM system. Due to our scaling of M, these costs actually decrease with T. 

5.2. Couvergeuce. The accuracy of the solution obtained through our algorithms is shown in 
Figures 5.2a through 5.2d. The accuracy is measured as the maximum pointwise error over an equi- 
sampled grid ten times denser than the one used for construction. This is measured for increasing 
number of degrees of freedom N for four test functions: 
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Fig. 5.1: Execution time for increasing degrees of freedom N, for MATLAB implementations of the 
Explicit and Implicit projection algorithms, the Lyon algorithm and a direct solver. 


• A well-behaved, smooth function to show convergence in near optimal conditions, 

fi{x) = x^. 

• A highly oscillatory function, to show the resolution power of Fourier extensions for oscillatory 
functions. 


f 2 {x) = Ai(76a;). 


• A function with a pole in the complex plane near the real interval [—1,1], to show convergence 
at a slower, but still exponential rate. 


h{x) 


1 

1.1-x"^' 


• A function with discontinuous first derivative, to see the breakdown of the algorithm to alge¬ 
braic convergence speeds, 


f 4 {x) = |a;|. 

The convergence behaviour seen is in accordance with [18], [2] and [3]. Convergence for functions 
analytic in [—1,1] is at least geometric, even when singularities are present near the real interval. 
Following the earlier arguments about resolution power from §2.3, Fourier extensions of oscillatory 
functions start to converge sooner for lower values of T. Note that this is in terms of degrees of 
freedom, and that we increased the oversampling for lower T to maintain conditioning. 

When the approximant is in C^, having continuous derivatives up to the convergence rate 
becomes algebraic at a rate of The convergence of order 0{N~^) for the (7° function is 

apparent from Figure 5.2d. 
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Fig. 5.2: The L^o norm of the error, computed by oversampling the solution by a factor 10, for the 
various testfunctions 


5.3. Robustness. To ensure there is no accumulation of error for very large IV, Figure 5.3 shows 
the continuation of the accuracy experiment for N up to 10^. No error accumulation is visible, the 
error fluctuates around the cutoff threshold. 

Figure 5.3a shows the challenging problem of approximating an oscillatory function with a fre¬ 
quency that increases with N, f{x) = sin{Nx). The error for the Algorithm using T = 1.1 and T = 2 
stays close to machine precision, if maybe slightly increasing. However, for T = 3.8, the maximum 
frequency in the Fourier basis is lower than the frequency of the signal for every N, so there cannot be 
any convergence. This experiment shows that even a very oscillatory signal such as f{x) = sin(10®a:) 
can be accurately approximated, as long as there are sufficient degrees of freedom for the chosen value 
of T. 
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